% jkl.m

function [PI,PHIs,PHIn, ...
    Ed2z,T2zq,T2z,T2qq,T2qv,T2q, ...
    Ed3z,T3zq,T3z,T3qq,T3qv,T3q] = jkl(varargin)

% Syntax
% ======= 
% 
% The code will provide a first-, second-, or third-order perturbation
% approximation solution based on numerical arrays of derivatives of the
% model equtaions. For example, 
%
%    -- first order:
%    [PI,PHIs,PHIn] = jkl(Gzp,Fzp,Gz,Fz,Gnp,Fnp, ...
%        Onn)
%
%    -- second order:
%    [PI,PHIs,PHIn, ...
%       Ed2z,T2zq,T2z,T2qq,T2qv,T2q] 
%       = jkl(Gzp,Fzp,Gz,Fz,Gnp,Fnp, ...
%	    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, ...
%       Onn)
%
%    -- third order:
%    [PI,PHIs,PHIn, ...
%       Ed2z,T2zq,T2z,T2qq,T2qv,T2q, ...
%       Ed3z,T3zq,T3z,T3qq,T3qv,T3q] 
%       = jkl(Gzp,Fzp,Gz,Fz,Gnp,Fnp, ...
%	    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, ...
%       GFzpzpzp,GFzpzpz,GFzpzpnp,GFzpzz,GFzpznp,GFzpnpnp,GFzzz,GFzznp,GFznpnp,GFnpnpnp,...
%       Onn)
% 
% Options
% ========
% -  'ignore=': ignore errors (default: false)

%% handle input arguments
[ignore,varargin] = parseOption(varargin,'ignore',false) ; % to ignore errors

it = numel(varargin) ;
if it == 7
    order = 1 ;
elseif it == 13
    order = 2 ;
elseif it == 23 ;
    order = 3 ;
else
    error('jkl: %s\n','Incorrect input argument syntax.') ;
end

Gzp = varargin{1} ;
Fzp = varargin{2} ;
Gz = varargin{3} ;
Fz = varargin{4} ;
Gnp = varargin{5} ;
Fnp = varargin{6} ;
if order>1
    GFzz = varargin{7} ;
    GFzpz = varargin{8} ;
    GFnpnp = varargin{9} ;
    GFznp = varargin{10} ;
    GFzpnp = varargin{11} ;
    GFzpzp = varargin{12} ;
    if order>2
        GFzpzpzp = varargin{13} ; %#ok<*NASGU>
        GFzpzpz = varargin{14} ;
        GFzpzpnp = varargin{15} ;
        GFzpzz = varargin{16} ;
        GFzpznp = varargin{17} ;
        GFzpnpnp = varargin{18} ;
        GFzzz = varargin{19} ;
        GFzznp = varargin{20} ;
        GFznpnp = varargin{21} ;
        GFnpnpnp = varargin{22} ;
    end
end
Onn = varargin{end} ;
nz = size(Fz,2) ;
nn = size(Fnp,2) ;

% check input arguments
if order>1
    assert(all(size(GFzz)==[nz,nz^2]),'Check input arguments.') ;
    assert(all(size(GFzpz)==[nz,nz^2]),'Check input arguments.') ;
    assert(all(size(GFnpnp)==[nz,nn^2]),'Check input arguments.') ;
    assert(all(size(GFznp)==[nz,nz*nn]),'Check input arguments.') ;
    assert(all(size(GFzpnp)==[nz,nz*nn]),'Check input arguments.') ;
    assert(all(size(GFzpzp)==[nz,nz^2]),'Check input arguments.') ;
    if order>2
        assert(all(size(GFzpzpzp)==[nz,nz^3]),'Check input arguments.') ;
        assert(all(size(GFzpzpz)==[nz,nz^3]),'Check input arguments.') ;
        assert(all(size(GFzpzpnp)==[nz,nz^2*nn]),'Check input arguments.') ;
        assert(all(size(GFzpzz)==[nz,nz^3]),'Check input arguments.') ;
        assert(all(size(GFzpznp)==[nz,nz^2*nn]),'Check input arguments.') ;
        assert(all(size(GFzpnpnp)==[nz,nz*nn^2]),'Check input arguments.') ;
        assert(all(size(GFzzz)==[nz,nz^3]),'Check input arguments.') ;
        assert(all(size(GFzznp)==[nz,nz^2*nn]),'Check input arguments.') ;
        assert(all(size(GFznpnp)==[nz,nz*nn^2]),'Check input arguments.') ;
        assert(all(size(GFnpnpnp)==[nz,nn^3]),'Check input arguments.') ;
    end
end

%% First order

% First order solution with dummy process
[exu,~,PHIs,~,~,PI,~,PHIn] = qzs(-Gzp,-Fzp,Gz,Fz,Gnp,Fnp,'quiet=',false) ;
if ~ignore
    if ~exu(1)
        error('jkl: %s\n','Solution does not exist.') ;
    end
    if ~exu(2)
        error('jkl: %s\n','Solution is not unique.') ;
    end
end
ns = size(PHIs,1) ;
Oss = reshape( (eye(ns^2)-kron(PHIs,PHIs))\kron(PHIn,PHIn)*vec(Onn), ns, ns ) ;
ng = size(Gz,1) ;
nf = size(Fz,1) ;
ny = ng + nf ;

%% Second order
if order>1
    % Dx is a duplication matrix such that Dx*vech(x*x') == vec(x*x')
    Dn=h2v(nn) ;
    Ds=h2v(ns) ;
    iDs = pinv(Ds) ;
    % Ksn is a commutation matrix
    Ksn = commute(ns, nn) ;
    
    nc2 = ns*(ns+1)/2 ;
    nv2 = nn*(nn+1)/2 + ns*nn ;
    nq2 = ns + ns*(ns+1)/2 ;
    
    %     C21 = [Fzpzp; Gzpzp]*kron(PI*PHIs,PI*PHIs)*Ds + 2*[Fzpz; Gzpz]*kron(PI*PHIs,PI)*Ds + [Fzz; Gzz]*kron(PI,PI)*Ds ;
    %     C22 = [Fzpzp; Gzpzp]*kron(PI*PHIn,PI*PHIn)*Dn + 2*[Fzpnp; Gzpnp]*kron(PI*PHIn,eye(nn))*Dn + [Fnpnp; Gnpnp]*Dn ;
    %     C23 = [Fzpzp; Gzpzp]*( kron(PI*PHIn,PI*PHIs) + kron(PI*PHIs,PI*PHIn)*Ksn ) ...
    %         + 2*[Fzpz; Gzpz]*kron(PI*PHIn,PI) + 2*[Fzpnp; Gzpnp]*kron(PI*PHIs,eye(nn))*Ksn ...
    %         + 2*[Fznp; Gznp]*kron(PI,eye(nn))*Ksn ;
    %     C2 = [C21, C22, C23] ;
    C2 = NaN(ny,nc2+nv2) ;
    C2(:,1:nc2) = GFzpzp*kron(PI*PHIs,PI*PHIs)*Ds + 2*GFzpz*kron(PI*PHIs,PI)*Ds + GFzz*kron(PI,PI)*Ds ;
    C2(:,nc2+1:nc2+nn*(nn+1)/2) = GFzpzp*kron(PI*PHIn,PI*PHIn)*Dn + 2*GFzpnp*kron(PI*PHIn,eye(nn))*Dn + GFnpnp*Dn ;
    C2(:,nc2+nn*(nn+1)/2+1:nc2+nv2) = GFzpzp*( kron(PI*PHIn,PI*PHIs) + kron(PI*PHIs,PI*PHIn)*Ksn ) ...
        + 2*GFzpz*kron(PI*PHIn,PI) + 2*GFzpnp*kron(PI*PHIs,eye(nn))*Ksn ...
        + 2*GFznp*kron(PI,eye(nn))*Ksn ;
    
    gam21 = [vech(Oss); ...
        vech(Onn); ...
        zeros(nn*ns,1)] ;
    gam22 = [eye(ns*(ns+1)/2)
        zeros(nn*(nn+1)/2+ns*nn,ns*(ns+1)/2)] ;
    gam23 = [zeros(ns*(ns+1)/2,nn*(nn+1)/2+nn*ns); ...
        eye(nn*(nn+1)/2) zeros(nn*(nn+1)/2,nn*ns); ...
        zeros(ns*nn,nn*(nn+1)/2) eye(ns*nn)] ;
    
    phi21 = zeros(nc2,1) ;
    phi22 = iDs*kron(PHIs,PHIs)*Ds ;
    phi23 = [iDs*kron(PHIn,PHIn)*Dn,iDs*(kron(PHIn,PHIs)+kron(PHIs,PHIn))*Ksn] ;
    
    % Second order solution
    % First use [exu,Tk,Tkk,Tkw,Ty,Tyk,Tyw,Tkv] = qzs() to compute fully
    % demeaned solution:
    [~,T2k,T2kk,T2kc,T2y,T2yk,T2yc,T2kv] = qzs(-Gzp,-Fzp,Gz,Fz,C2(1:ng,:),C2(ng+1:ny,:),...
        gam21,gam22,gam23,phi21,phi22,phi23,'useoldqz=',true) ;
    
    % Then compute non-demeaned solution:
    Ed2s = (eye(ns)-PHIs)\T2k ; % unconditional expectation of second order differential
    Ed2z = T2y + T2yk*Ed2s ;
    
    T2zq = [T2yk T2yc] ;
    T2z = Ed2z - T2zq*[Ed2s; vech(Oss)] ;
    T2qq = [T2kk, T2kc; zeros(nc2,ns), phi22] ;
    T2qv = [T2kv; phi23] ;
    T2q = (eye(size(T2qq))-T2qq)*[Ed2s; vech(Oss)] - T2qv*[vech(Onn); zeros(ns*nn,1)] ;
    
    
    %% Third order
    if order>2
        %         clear C2 gam2* phi2*
        Kzz = commute(nz,nz) ;
        
        nc3 = ns + ns*nq2 ;
        
        g_zpzp = [zeros(nz^2,ns) kron(PI*PHIs,PI*PHIs)*Ds] ;
        a_zpzp = [kron(PI*PHIn,PI*PHIn)*Dn, kron(PI*PHIn,PI*PHIs)+kron(PI*PHIs,PI*PHIn)] ;
        g_zpz = [zeros(nz^2,ns), kron(PI*PHIs,PI)*Ds] ;
        a_zpz = [zeros(nz^2,nn*(nn+1)/2), kron(PI*PHIn,PI)] ;
        a_zpnp = [kron(PI*PHIn,eye(nn))*Dn kron(PI*PHIs,eye(nn))*Ksn] ;
        g_zz = [zeros(nz^2,ns), kron(PI,PI)*Ds] ;
        a_znp = [zeros(nz*nn,nn*(nn+1)/2), kron(PI,eye(nn))*Ksn] ;
        a_npnp = [Dn, zeros(nn^2,nn*ns)] ;
        
        %         C3t = NaN(ny,ns+nq2+nn+nn*nq2+ns*nv2+nn*nv2) ;
        C31 = 3*GFzpzp*( kron(T2z,PI)*PHIs + kron(T2zq*T2q,PI*PHIs) ) ...
            + 3*GFzpz*( kron(T2z,PI) + kron(T2zq*T2q,PI) + Kzz*kron(T2z,PI*PHIs) ) ...
            + 3*GFzz*kron(T2z,PI) ;
        C32 = 3*GFzpzp*kron(T2zq*T2qq,PI*PHIs) ...
            + 3*GFzpz*( kron(T2zq*T2qq,PI) + Kzz*kron(T2zq,PI*PHIs) ) ...
            + 3*GFzz*kron(T2zq,PI) ...
            + GFzpzpzp*kron(g_zpzp,PI*PHIs) ...
            + GFzpzpz*kron(g_zpzp,PI) ...
            + GFzpzz*kron(g_zpz,PI) ...
            + GFzzz*kron(g_zz,PI) ;
        C33 = 3*GFzpzp*( kron(T2z,PI)*PHIn + kron(T2zq*T2q,PI*PHIn) ) ...
            + 3*GFzpnp*( kron(T2z,eye(nn)) + kron(T2zq*T2q,eye(nn)) ) ...
            + 3*GFzpz*Kzz*kron(T2z,PI*PHIn) ...
            + 3*GFznp*kron(T2z,eye(nn)) ;
        C34 = 3*GFzpzp*kron(T2zq*T2qq,PI*PHIn) ...
            + 3*GFzpnp*kron(T2zq*T2qq,eye(nn)) ...
            + 3*GFzpz*Kzz*kron(T2zq,PI*PHIn) ...
            + 3*GFznp*kron(T2zq,eye(nn)) ...
            + GFzpzpzp*kron(g_zpzp,PI*PHIn) ...
            + 3*GFzpzpnp*kron(g_zpzp,eye(nn)) ...
            + 6*GFzpznp*kron(g_zpz,eye(nn)) ...
            + 3*GFzznp*kron(g_zz,eye(nn)) ;
        C35 = 3*GFzpzp*kron(T2zq*T2qv,PI*PHIs) ...
            + 3*GFzpz*kron(T2zq*T2qv,PI) ...
            + GFzpzpzp*kron(a_zpzp,PI*PHIs) ...
            + 3*GFzpzpz*kron(a_zpzp,PI) ...
            + 3*GFzpzz*kron(a_zpz,PI) ;
        C36 = 3*GFzpzp*kron(T2zq*T2qv,PI*PHIn) ...
            + 3*GFzpnp*kron(T2zq*T2qv,eye(nn)) ...
            + GFzpzpzp*kron(a_zpzp,PI*PHIn) ...
            + 3*GFzpzpnp*kron(a_zpzp,eye(nn)) ...
            + 6*GFzpznp*kron(a_zpz,eye(nn)) ...
            + 3*GFzpnpnp*kron(a_zpnp,eye(nn)) ...
            + 3*GFznpnp*kron(a_znp,eye(nn)) ...
            + GFnpnpnp*kron(a_npnp,eye(nn)) ;
        C3 = [C31, C32, C33, C34, C35, C36] ;
        
        gam31 = [zeros(ns,1); ...
            zeros(ns*nq2,1); ...
            zeros(nn+nn*nq2+ns*nv2+nn*nv2,1)] ;
        gam32 = [eye(ns), zeros(ns,ns*nq2); ...
            zeros(ns*nq2,ns), eye(ns*nq2); ...
            zeros(nn+nn*nq2+ns*nv2+nn*nv2,ns+ns*nq2)] ;
        gam33 = [zeros(ns+ns*nq2,nn+nn*nq2+ns*nv2+nn*nv2); ...
            eye(nn), zeros(nn,nn*nq2+ns*nv2+nn*nv2); ...
            zeros(nn*nq2,nn), eye(nn*nq2), zeros(nn*nq2,ns*nv2+nn*nv2); ...
            zeros(ns*nv2,nn+nn*nq2), eye(ns*nv2), zeros(ns*nv2,nn*nv2); ...
            zeros(nn*nv2,nn+nn*nq2+ns*nv2), eye(nn*nv2)] ;
        
        nc3 = ns + ns*nq2 ;
        nv3 = nn + nn*nq2 + ns*nv2 + nn*nv2 ;
        
        phi31 = zeros(nc3,1) ;
        phi32 = [PHIs zeros(ns,ns*nq2); ...
            kron(T2q,PHIs) kron(T2qq,PHIs)] ;
        phi33 = [PHIn zeros(ns,nn*nq2+ns*nv2+nn*nv2); ...
            kron(T2q,PHIn) kron(T2qq,PHIn) kron(T2qv,PHIs) kron(T2qv,PHIn)] ;
        
        % First use [exu,Tk,Tkk,Tkw,Ty,Tyk,Tyw,Tkv] = qzs() to compute fully
        % demeaned solution:
        [~,T3k,T3kk,T3kc,T3y,T3yk,T3yc,T3kv] = qzs(-Gzp,-Fzp,Gz,Fz,C3(1:ng,:),C3(ng+1:ny,:),...
            gam31,gam32,gam33,phi31,phi32,phi33,'useoldqz=',true) ;
        
        % Third order solution
        Ed3s = (eye(ns)-T3kk)\T3k ; % unconditional expectation of second order differential
        Ed3z = T3y + T3yk*Ed3s ;
        
        T3zq = [T3yk T3yc] ;
        T3z = Ed3z - T3zq*[Ed3s; zeros(nc3,1)] ;
        T3qq = [T3kk, T3kc; zeros(nc3,ns), phi32] ;
        T3qv = [T3kv; phi33] ;
        T3q = (eye(size(T3qq))-T3qq)*[Ed3s; zeros(nc3,1)] ;
        
    end % if order>2
end % if order>1
end % function 





